REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No.  0704-0188 

The  public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources, 

gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of 
information,  including  suggestions  for  reducing  the  burden,  to  the  Department  of  Defense,  Executive  Services  and  Communications  Directorate  (0704-0188).  Respondents  should  be  aware 
that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection  of  information  if  it  does  not  display  a  currently  valid  OMB 
control  number. 

PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ORGANIZATION. 

1.  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE 

20  Mar  2008  FINAL  REPORT 

3.  DATES  COVERED  (From  -  To) 

1  Mar  2004  -  3 1  Jul  2007 

4.  TITLE  AND  SUBTITLE 

QUANTUM  LATTICE  ALGORITHMS  FOR  2D  AND  3D 
MAGNETOHYDRODYNAMICS 

5a.  CONTRACT  NUMBER 

FA9550-04- 1-0128 

5b.  GRANT  NUMBER 

03NM208 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

Dr.  Linda  Vahala 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Dept  of  Electrical  and  Computer  Engineering 

OLD  DOMINION  UNIV  RSH  FOUNDATION  INC. 

Norfolk,  VA  23508 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

AF  OFFICE  OF  SCIENTIFIC  RESEARCH 

875  NORTH  RANDOLPH  STREET  ROOM  31 12 

ARLINGTON  VA  22203  . 

10.  SPONSOR/MONITOR  S  ACRONYM(S) 

11.  SPONSOR/MONITOR  S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

DISTRIBUTION  STATEMENT  A:  UNLIMITED 


AFRL-SR-AR-TR-08-01 76 


13.  SUPPLEMENTARY  NOTES 


14.  ABSTRACT 

During  the  3  years  of  this  grant,  we  have  continued  our  collaboration  with  Jeffrey  Yepez 
(AFRL,  Hancom  Field)  and  George  Vahala  (William  &  Mary)  on  both  quantum  and 
entropic  lattice  algorithms  for  the  solution  of  nonlinear  physics  problems.  Because  of  the 
extreme  scalability  of  the  algorithms  that  we  have  been  developing,  we  were  chosen  for 
CAP-Phase  II  for  the  new  IBM-P5+  supercomputer  (Babbage)  at  NAVO  MSRC  and  also 
for  CAP-Phase  II  on  the  9000  core  on  the  SGI-Altix  at  ASC. 


15.  SUBJECT  TERMS 

Nonlinear  Physics;  Quantum  Lattice  Algorithms;  Entropic  Lattice  Algorithms 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

ABSTRACT 

18.  NUMBER 

OF 

PAGES 

14 

19a.  NAME  OF  RESPONSIBLE  PERSON 

a.  REPORT 

b.  ABSTRACT 

c.  THIS  PAGE 

19b.  TELEPHONE  NUMBER  (Include  area  code) 

Standard  Form  298  (Rev.  8/98) 

Prescribed  by  ANSI  Std.  Z39.18 


AFOSR  FINAL  REPORT 


"QUANTUM  LATTICE  ALGORITHMS  FOR  2D  AND  3D 
MAGNETOHYDRODYNAMICS" 


November  2007 


Linda  Vahala 
Old  Dominion  University 


200804041 20 


During  the  3  years  of  this  grant,  we  have  continued  our  collaboration  with  Jeffrey  Yepez 
(AFRL,  Hancom  Field)  and  George  Vahala  (William  &  Mary)  on  both  quantum  and 
entropic  lattice  algorithms  for  the  solution  of  nonlinear  physics  problems.  Because  of  the 
extreme  scalability  of  the  algorithms  that  we  have  been  developing,  we  were  chosen  for 
CAP-Phase  II  for  the  new  IBM-P5+  supercomputer  (Babbage)  at  NAVO  MSRC  and  also 
for  CAP-Phase  II  on  the  9000  core  on  the  SGI-Altix  at  ASC. 


What  is  very  interesting  is  the  analogy  between  the  detailed  balance  quantum 
lattice  algorithms  and  entropic  lattice  Boltzmann  algorithms. 

At  each  space-time  grid  point  (x,f)  in  lattice  algorithms,  the  excited  state  of  a  qubit  |  <7^ 
encodes  the  probability  of  the  existence  of  a  mesoparticle  moving  with  discrete  lattice 
velocity  c?  =  Ax^  /  At .  Ax^are  the  lattice  vector  links,  with  q  =  1,2, ....  Q ,  where  Q  is 

the  number  of  qubits  at  each  spatial  node.  The  particle  momentum  is  determined  from  a 
suitably  chosen  qubit-qubit  interaction  Hamiltonian  H'  while  the  spatial  location  arises 
from  the  free-streaming  Hamiltonian  -ih£ V.  All  the  particle-particle  interactions 

9 

generated  by  H'  (from  2-body  up  to  Q-body  interactions)  can  be  mapped  onto  a  local 
collision  operator  | j  at  x.  In  particular,  for  type-II  quantum  algorithms,  the 

quantum  entanglement  is  localized  to  those  Q-qubits  at  (x,/)  and  then  this  entanglement 
is  spread  throughout  the  lattice  by  unitary  streaming3,4: 

=  +  >  4(x+av+a')=4'M  •  o) 


Here  is  the  incoming  probability  and  the  outgoing  probability.  In  the  classical 
limit,  there  exists  a  fundamental  discrete  entropy  function1 ,2,5 

?=> 


where  the  normalized  weights 


\ 

y  w  =1 
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are  determined  self-consistently. 


The 


9  J 

collision  operator  in  Eq.  (1)  is  determined  so  that  one  remains  on  a  constant  entropy 


surface 

Hi.','  -  « (/.../„)  •  (3) 

Eqs.  (l)-(3)  constitute  the  basics  of  the  detailed-balance  lattice  algorithms  for  fluid 
turbulence  that  are  ideal  for  parallel  (both  classical  and  quantum)  supercomputers. 

In  the  Q-dimensional  velocity  space,  the  relaxation  distribution  function  feq  is 
determined  analytically  by  extremizing  the  H-function  subject  to  the  local  collisional 
constraints  of  conservation  of  probability  and  probability  flux.  f*q ,  considered  as  a 

vector,  is  the  bisector  of  the  difference  between  the  incoming  and  outgoing  kinetic 
vectors  in  the  inviscid  limit  lim/(^0  a  /  2r  =  2  : 


(2) 


Eliminating  Q^and  fq 'from  Eqs.  (4)  and  (1)  one  obtains  the  lattice  Boltzmann  (LB) 
equation 


/,(x  +  Ai,,(  +  A/)=/t(x,<)  +  ^[/^(x,/)- /,(!,«)]  ,  q  =  \....Q  (5) 

This  is  basically  the  entropic  LB1,2  with  the  BGK  collisional  relaxation  parameters 
a(x,/)/2rand  feq  determined  from  Eqs.  (2)  and  (3).  In  the  Chapman-Enksog  limit, 

(Ax-h>0  ,  At  — >  o)—  and  identifying  the  density  and  momentum  moments  =  p, 
£c  /  =  p  u  —  one  recovers  the  quasi-incompressible  Navier-Stokes  equation  with 


effective  viscosity:  p(x,r)  =  — 


4t 


a(x,f) 


-  1 


molecular  viscosity:  p0  = -(2r- l)  ,  t>0.5 


(6) 


To  avoid  discrete  lattice  geometry  effects  polluting  the  turbulence  simulations,  one  is 
restricted  to  certain  Q’s  on  a  cubic  lattice.  In  particular  it  can  be  shown  that  on  a  unit 
cubic  lattice,  the  lowest  order  kinetic  velocity  models  are 

Q15:  rest  velocity,  speed  1  (6  velocities),  speed  Vi  (8  velocities)  -  i.e.,  Q=  15 

Q19:  rest  velocity,  speed  1  (6  velocities),  speed  Vi  (12  velocities)  -  i.e.,  0=19 

Q27:  rest  velocity,  speed  1  (6),  speed  Vi  (12),  and  speed  Vi  (8)  -  i.e.,  Q  =  27  (7) 

Because  detailed  balance  is  in-built  into  the  entropic  LB  algorithm  [see  Eq.  (3)],  the 
scheme  is  unconditionally  stable  for  arbitrary  large  Reynolds  numbers.  Re  =  UQL  /  2/rp0 . 


In  the  CAP-II  runs  on  9000  cores  on  Hawk,  we  concentrated  on  the  Lattice  Boltzman 
MHD  code.  In  particular  we  have  been  considering  the  lattice  Boltzmann  simulation  for 
the  initial  profiles  of  a  Taylor-Green  velocity  profile  in  an  Orszag-Tang  magnetic  field: 
u(x,  t  =  0)  =  U0  (sinxcosycosz  ,  -cosxsinycosz,  0) 

B(x,  t  =  0)  =  B0  (-2sin2y  +  sinz,  2sinx  + sinz,  sinx  +  siny) 

The  isosurfaces  of  vorticity  and  current  shown  in  Fig.  1 


Fig.  1  (a)  Isosurfaces  of  \Vorticity\  (b)  Isosurfaces  of  \Current\ 

With  these  profiles,  there  is  no  initial  magnetic  helicity  or  cross  helicity 

0  =  J d3* A(x,0)»B(x,0)  ,  0  =  Jd3;tu(x,0)»B(x,0) 

where  A  is  the  vector  potential. 

An  extremely  important  property  of  the  LB  algorithm  for  MHD  is  that  from  the 
Chapman-Enksog  expansions  one  can  show  that  the  trace  of  the  first  order  magnetic 
stress  tensor  is  proportional  the  divergence  of  the  magnetic  field  -  and  hence  this  must  be 
zero  since  the  magnetic  stress  tensor  is  antisymmetric: 

0  =  7>A(1)  =  2>«,  [gai  -C]  =  -^V*B 

We  have  verified  this  result  directly  from  our  LB  simulations  by  explicitly  calculating  the 
trace  of  the  magnetic  stress  tensor,  Fig.  2,  with  7VA* ' 1  =  0  to  machine  accuracy. 


Maxt  !  .V  -  B  ( * ,  / )  I  /  B0 


Fig. 2  The  time  evolution  of  TrA(t)(t)  =  £eai(gai  -  g%)  in  the  LB  simulation,  showing  it  is  0  to 

a,  i 

machine  accuracy.  Chapman-Enskog  asymptotics  yields  TrA(l)  =  -TbV*B  /  3 . 

We  present  results  from  two  large  simulation  runs  on  a  spatial  lattice  of 
1 800  x  1 800  x  1 800  using  all  the  9000  cores  available  on  the  SGI  Altix.  The  first  run  ran 
for  60  K  time  steps  with 

CASE  A  :  Re  =  ^l  =  1000  ,  Rm  =  ^  =  350  ,  Pr  =  -  =  0.3  , 
v  r\  ii 


while  the  second  run  ran  to  30K  time  steps  at  a  higher  Prandtl  and  magnetic  Reynolds 
number 

CASE  B  :  Re  =  =  350  ,  Rm  =  ^-  =  1050,  Pr  =  -=3.0 

v  t?  n 

In  Fig.  3  we  plot  the  time  development  of  the  normalized  energies,  enstrophies  and 
palinstrophy 

£*ta(0=Jrf3*u2M  ’  0=K*B2M 

Kinetic  Enstrophy  Cl(t)  =  J d*x  |V  x  u(x,f)|  =  ((02(x,t)^ 

Magnetic  Enstrophy  QM(t)  =  J d3x  |V  xB(x,/)|  =  ^J“(x,t)^ 

Palinstrophy  P(t)  =  jd*x  |V X(o(x,r)|  =  ^|V XQ)(x,r)|  ^ 

which  are  just  higher  order  k—  moments  of  the  energy  spectra. 

Normalized  Enstropy 


Normalized  Palinstrophy 


o 
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(c)  normalized  palinstrophy 

Fig.  3  The  time  development  of  the  (a)  kinetic  and  magnetic  energies,  (b)  the  kinetic  and 
magnetic  enstrophy,  and  (c)  the  kinetic  palinstrophy  -  normalized  to  peak  value. 


In  Navier-Stokes  turbulence,  the  kinetic  enstrophy  increases  sharply  at  early  times  due  to 
inviscid  vortex  stretching  -  i.e.,  the  kinetic  enstrophy  increase  is  independent  of  the 
transport  coefficient.  The  kinetic  energy  during  this  period  is  very  slowly  decreasing.  In 
MHD,  however,  we  see  an  immediate  strong  energy  exchange  from  kinetic  to  magnetic 
which  is  independent  of  transport  coefficients  (Case  A  and  B  curves  overlay  in  Fig.  4a  for 
t  <  10  K)  with  a  rapid  rise  in  the  magnetic  enstrophy  (i.e.,  mean  square  current).  For  10K 
<  t  <  20K,  there  is  a  flattening  in  the  kinetic  energy  decay  (Fig.  4a)  and  a  subsequent 
increase  in  the  kinetic  enstrophy  (Fig.  4b),  somewhat  akin  to  Navier-Stokes  turbulence. 
The  strength  of  the  respective  transport  coefficient  dictates  which  particular  enstrophy 
peaks  at  a  greater  value  (i.e.,  for  Case  B  the  lower  resistivity  and  higher  viscosity  dictate 
that  the  magnetic  enstrophy  has  a  greater  increase  than  in  Case  A  while  the  kinetic 
enstrophy  has  a  lower  increase  than  in  Case  A).  This  is  also  seen  in  the  sharp  rise  of  the 
kinetic  palinstrophy,  Fig.  4c. 

The  directional  energy  spectra  are  shown  in  Fig.  5  (low  magnetic  Prandtl  number, 
Pr  =  v/77  =  0.3)  and  Fig.6  (high  magnetic  Prandtl  number,  Pr  =  v / rj  =  3.0 ).  Initially 

these  spectra  are  delta  functions.  The  directional  kinetic  and  magnetic  spectra  are  defined 
by 


2 

^  \llx ,  EKy(kx,t)  —  ^  \uy (kx,ky,kZ',t^ 


k,  k 


ky ,  k.  ky ,  k. 

where  the  summation  is  always  over  the  wavenumbers  ky  ,  k, :  the  longitudinal  spectra 
involve  the  x-component  of  the  fields  while  the  transverse  spectra  the  y-component  of  the 
fields.  In  Fig.  4a,  we  plot  the  longitudinal  [initially  a  <5(^-2)-  spectrum]  and  the 
transverse  [initially  a  two-delta  function  peak  spectrum  at  kx=  2,  4  ]  directional  kinetic 
energy  at  time  =  30K,  while  in  Fig.  4b  the  directional  longitudinal  magnetic  energy 
spectrum  at  t  =  10K,  20K,  and  30K  -  and  the  comparison  to  the  k -;>/3  Kolmogorov 
spectrum. 


Fig.  4  The  time  development  of  the  directional  energy  spectra  for  the  low  magnetic  Prandtl 
number  case:  pr  =  v  /  r\  =  0.3  (a)  the  longitudinal  and  transverse  kinetic  energy  spectrum  at  t  = 

30K,  and  (b)  the  longitudinal  magnetic  energy  spectrum  at  t  =  10K,  20K  and  30K.  Also  plotted  is 
the  Kolmogorov  k'5/}  inertial  range  spectrum. 


number  case:  pr  =  v  /  fj  =  3.0  (a)  the  longitudinal  and  transverse  kinetic  energy  spectrum  at  t  - 
4K  and  t  -  16K,  and  (b)  the  longitudinal  and  transverse  magnetic  kinetic  energy  spectrum  at  t  = 
4K  and  t  =  16K.  Also  plotted  is  the  Kolmogorov  k'5/}  inertial  range  spectrum. 

As  we  increase  the  magnetic  Prandtl  number  to  Case  B  one  finds  a  substantial  difference 
between  the  longitudinal  and  transverse  spectra.  Fig.  5.  The  transverse  kinetic  and 
magnetic  energy  spectra  show  much  stronger  excitation  of  high  kx  -modes,  which  in  the 
magnetic  energy  case  shows  a  quite  strong  semblance  to  the  Kolmogorov  inertial  energy 


k~ 5/3  -  spectrum  at  t  =  16K  (Fig.  5b).  There  is  also  an  interesting  enhancement  of  the 
transverse  kinetic  energy  spectrum  for  40  <  kx  <  200  at  t  =  16K  (Fig.  5a). 

This  behavior  maybe  attributed  to  the  strong  but  localizd  vorticity  and  current 
sheets  developing  due  to  the  turbulence.  In  Fig.  6  we  plot  some  time  snapshots  of  the 
isosurfaces  of  vorticity  and  current  for  this  high  magnetic  Prandtl  number  case. 


vorticity  isosurface 


current  isosurface 


current  isosurface 


vorticity  isosurface 


(a)  t  =  4K 


(b)  t  =  8K 


vorticity  isosurface 


vorticity  isosurface 


current  isosurface 


current  isosurface 


(c)  t  =  12K 


(d)  t  =  16K 


current  isosurface  current  isosurface 

(e)  t  =  20K  (f)  t=  24K 

Fig.  6  The  time  development  of  the  vorticty,  |d)| ,  and  corresponding  current  |  j| ,  isosurface  for 
the  high  magnetic  Prandtl  number  simulation,  (a)  t  =  4K,  (b)  t  =  8K,  (c)  t  —  I2K,  (d)  t  —  16K,  (e) 
t  =  20K,  (f)  t  =  24K.  The  isosurface  value  chosen  is  that  corresponding  to  the  average  |co|  and 

A  A  A 

average  |  j|  for  that  time  instant.  The  color  coding  is  dependent  on  the  value  of  U*CD  and  B*J  at 
the  isosurface  gridpoint,  going  from  RED  for  parallel  unit  vectors  U  ||  <0  to  BLUE  for 
antiparallel  unit  vectors  u*©=  -1.  Similarly  for  the  current  isosurfaces:  from  RED  for 
B*J  =  +1  to  BLUE  for  B»j  =  -1 .  The  GREY  scale  is  for  isosurfaces  with  U*G)  =  0  =  B*J 

There  is  much  information  in  Fig.  6  :  the  intensification  of  localized  horizontal  current 
sheets  (see  Fig.  6a-Fig.6c,  midway  at  the  vertical  cube  edges),  the  development  of  intense 
vertical  localized  patches  of  vorticity  and  current  at  later  times  with  similar  isosurface 


vorticity  isosurface 


vorticity  isosurface 


geometrical  structures  of  vorticity  and  current.  It  is  also  very  apparent  that  large  scale 
magnetic  (and  hence  velocity)  structures  persist  for  long  times,  Fig.6f.  This  is  also  seen 
in  the  low  magnetic  Prandtl  number  simulation,  Fig.  7,  where  some  large  scale  vorticity 
and  current  isosurface  structures  persist,  even  at  t  =  50K.  This  is  very  unlike  3D  Navier- 
Stokes  turbulence  which  is  dominated  by  small  scale  vortex  structures  as  seen  in  Fig.  8 


(a)  (b) 

Fig.  7  The  (a)  vorticity  and  (b)  current  isosurfaces  at  t  =  5 OK  for  the  low  Prandtl  number 
(Pr  =  v/rj  =  0.3)  simulation.  Some  large  scale  magnetic  structures  persist,  along  with 

corresponding  large  scale  vortex  structures. 


Fig.  8  The  vorticity  isosurfaces  in  3D  Navier-Stokes  turbulence.  The  flow  is  dominated  by  very 
small  scale  structures  after  the  inviscid  vortex  stretching  and  the  peak  in  the  fluid  enstrophy. 
This  isosurface  is  at  t  =  7 K  of  an  ELB  simulation 

Finally,  we  present  some  correlation  data  for  the  magnetic  field 
C,„ng(r)  =  {Bx(x,y,z)  Bx(x  +  r,y,z)) 


CLns(r)=(By(x’y’Z)  By{x  +  r,y,z))  ,  C?rans(r)  =  (Bz(x,y,z)  Bz{x  +  r,y,z)) 


Magnetic  Correlations  at  t  =  1 2  K  Magnetic  Correlations  at  t  =  24  K 


Fig.  14  Magnetic  Correlations  at  (a)  t  =  12K,  (b)  t  =  24K  for  the  high  magnetic  Prandtl 
number  simulation,  Case  B.  The  very  slight  increase  in  the  longitudinal  correlation  function  for  r 
>  650  at  t  =  12K  is  no  longer  present  at  t  =  24K. 

It  is  seen  that  as  time  develops,  the  almost  constant  asymptotic  tail  of  the  longitudinal 
magnetic  correlation  function  disappears,  and  by  t  =  28K  Clong  (r)  <  0  f°r  all  r 
Moreover,  we  find  (for  r  >  0 ) 

CLM  ,  C^ans{r)  <  Clong(r)  ,  for  all  times 

with  Clon'(r)  <  0,  are  consistent  with  the  correlation  statistics  of  a  random  solenoidal 
vector  field . 


Concluding  Remarks 

We  have  developed  a  3D  LB-MHD  algorithm  that  is  ideally  parallelized  and 
presented  some  simulation  results  on  a  1 8003-spatial  grid  that  shows  the  persistence  of 
large  scale  magnetic  and  vorticity  structures  for  long  times.  Morever,  the  time 
development  of  the  correlation  statistics  of  the  magnetic  field  indicate  that  the  B-field  is 
becoming  more  and  more  random.  An  important  feature  of  the  LB-MHD  approach  is 
that  the  algorithm  automatically  ensures  V*B  =  0  to  machine  accuracy. 

The  straightforward  LB  algorithm,  while  simple  and  explicit,  suffers  from 
numerical  instabilities  as  Re  ->  °°  ,  Rm  °°  .  This  places  upper  bounds  on  the 
attainable  transport  coefficients.  At  the  Navier-Stokes  level,  entropic  algorithms  have 
been  developed  that  remain  unconditionally  stable  for  arbitrary  small  viscosities. 
Indeed,  we  have  presented  here  the  first  large  scale  ELB  simulations  on  a  16003-grid  at 
Re  =  25000.  While  our  ELB  code  runs  successfully  for  much  higher  Reynolds  numbers, 
the  turbulence  is  no  longer  fully  resolved  on  these  ‘small’  grids,  and  so  these  results  are 
not  presented  here.  We  are  currently  developing  entropic  LB-MHD  algorithms  that 
would  permit  simulations  at  arbitrary  small  transport  coefficients. 


While  the  simulations  reported  here  are  on  a  simple  3D  periodic  domain,  LB 
algorithms  can  handle  arbitrary  geometries  without  loosing  their  intrinsic  parallelization. 
Nonuniform  spatial  grids  can  be  readily  handled.  In  these  cases,  the  spatial  grid  and  the 
kinetic  velocity  lattice  will  now  no  longer  overlay.  As  a  result,  the  streaming  step  of  the 
LB  algorithm  will  no  longer  give  immediate  data  at  the  spatial  nodes.  One  would  then 
resort  to  interpolation  methods  to  get  the  streamed  information  onto  the  spatial  nodes. 
Moreover  all  the  latest  CFD  methods  for  handling  arbitrary  spatial  grid  geometries  can  be 
immediately  brought  over  to  LB.  It  remains  to  be  seen  what  price  will  need  to  be  paid  on 
the  parallelization  of  such  augmented  LB  codes. 

Finally,  we  comment  on  another  interesting  aspect  of  ELB  algorithms.  The  ELB- 
viscosity  veff(\,t)  gives  the  appearance  of  an  eddy  viscosity  and  immediately  raises  the 

question  of  whether  there  is  any  connection  between  ELB  and  the  Large  Eddy 
Simulations  (LES)  in  turbulence  modeling.  The  simplest  LES  model  is  the  Smagorinsky 
model  in  which  the  subgrid  scales  are  modeled  by  an  eddy  viscosity  that  is  related  to  the 
mean  rate  of  strain  velocity  tensor: 

=  (C,A)2 


where  the  rate  of  strain  tensor  (of  the  resolvable  scales) 

_  1  3m,  ^  3 Uj 

"  2  V  dxj  +  3x,  J 

A  is  the  filter  width  (defined  in  the  filtering  function  that  separates  the  resolvable  from 
the  subgrid  scales)  and  Cs  is  some  empirical  constant.  Obviously,  the  connection  (if 
any)  between  the  ELB  and  LES  transport  coefficients  is  not  obvious:  ELB  deals  with 
entropy  surfaces  and  the  determination  of  the  collision  parameter  y(\,t)  that  enforces 

detailed  balance  on  the  pre-  and  post-collision  distribution  functions,  while  LES  deals 
with  the  rate  of  strain  tensor.  It  is  of  much  interest  that  one  can  immediately  construct 
local  LB-LES  models  that  recover  the  Smagorinsky-CFD  LES  model.  This  is  because 
the  local  strain  tensor  can  be  recovered  from  the  second  moment  of  the  non-equilibrium 
distribution  function 
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Of  course,  this  is  exactly  how  V*B  is  recovered  from  the  trace  of  the  first  moment  of  the 
nonequilibrium  magnetic  distribution  function  and  by  making  this  first  moment 

antisymmetric  we  enforce  V*B  =  0  to  better  than  O ( 1 0” 1 5 ) .  This  also  opens  up  the 


possibility  of  examining  LES  LB-MHD  algorithms  being  developed  for  CFD  techniques 
by  Carati  et.  al.,  where  the  LB  version  will  be,  unlike  the  CFD  code,  ideally  parallelized. 
These  LES  LB-MHD  codes  are  currently  being  developed. 
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